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Abstract 


A method  is  pro|io«v*d  for  determining  the  motion  of 
a I tody  relative  to  a fixrd  er\„onmcnl  using  the  changing 
image  seen  by  n inmcra  altc'lictl  to  the  l wily  The  optical 
flow  in  the  image  plane  i<  lb.  input,  while  the  instantaneous 
rotation  ami  translation  of  the  Ixxly  are  the  output.  If  op* 
tie  .1  How  could  lie  determined  precisely,  it  would  only  have 
to  l>e  kmiwu  at  a few  places  to  compute  the  parameter!  of 
the  motion.  In  practice,  however,  the  measured  optical  flow 
will  lie  wine wlmt  inaccurate  It  is  therefore  advantageous 
to  consider  methods  which  use  «ji  much  of  the  available  in* 
formation  a*  ]>o*>ible.  We  employ  a least* square*  approach 
which  miiiimues  some  measure  of  the  discrepancy  between 
tbe  measured  flow  and  that  predicted  from  the  computed 
motion  parameters  Several  different  error  norms  are  inves- 
tigated. In  general,  our  algorithm  leads  to  a system  of  non- 
liucar  equations  from  which  (he  motion  parameters  may  he 
computed  numerically.  However,  in  the  special  cases  where 
the  motion  of  the  camera  is  purely  iranrlaliotial  or  purely 
rotational,  uk  of  the  appropriate  norm  leads  to  a system  of 
equations  from  which  these  parameters  can  lie  determined 
in  closed  form,  v*- 
\ 
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1.  Introduction 

In  this  paper  we  investigate  the  problem  of  passive 
navigation  usiug  optical  flow  information.  Suppose  we  are 
viewing  n Dim.  We  wish  to  determine  the  motion  of  the 
camera  from  the  sequence  of  images,  assuming  that  the  in- 
stantaneous velocity  of  the  brightness  patterns,  also  called 


the  optical  How,  is  known  at  each  point  iu  the  image. 
Several  Kbemata  for  computing  optical  Dow  have  been  sug- 
gested (e.g.  |2|,  [3],  (5|).  Other  papers  (e.g.  (9),  (I  I],  (12)) 
have  previously  addresser!  the  problem  of  passive  naviga- 
tion. Three  approaches  can  be  taken  towards  a solution 
which  wc  term  the  discrete,  the  differentia/  and  the  con- 
tinuous approach. 

In  the  discrete  approach,  information  about  the  move- 
ment of  brightness  pattern*  at  only  a few  point*  is  used  to 
determine  the  motion  of  the  camera.  In  particular,  using 
such  an  approach,  one  attempts  lo  identify  and  match  dis- 
crete |H)inU  in  a sequence  of  images  Of  interest  iu  this 
case  is  the  photogramiiictric  problem  of  determining  what 
tbe  minimum  number  of  {mints  is  from  which  the  motion 
can  be  calculated  for  a given  number  of  images -(10,  (ll), 
(12),  JIG),  (17).  This  approach  requires  that  ouc  tracks  fea- 
tures, or  identifies  corresponding  feature*  in  images  taken 
at  different  times.  In  their  work,  Tsai  and  Ilunag  (IB)  as- 
sumed that  such  corresponding  {mints  ran  be  determined  in 
two  image.  Then  they  showed  that  in  grueral  seven  points 
are  suilicicut  to  determine  the  motion  uniquely.  They  prove 
furthermore  that  such  {mints  have  to  satisfy  a fairly  weak 
constraint.  Longuet-lliggins  work  (10)  is  fairly  similar  to 
(1C)  but  lie  fails  to  show  under  which  conditions  the  motion 
can  be  determined  uniquely  from  corresponding  points. 

In  tbe  differential  approach,  the  lire l and  Kcond  spatial 
partial  derivatives  of  the  optical  flow  are  used  to  compute 
tbe  motiou  of  a camera  [6],  (9).  It  has  been  claimed  that  it 
is  sufficient  to  know  the  optical  flow  and  both  its  first  and 
second  derivatives  at  a single  point  to  uniquely  determine 
the  motion  (9).  This  is  incorrect  (except  for  a special  case) 
(1).  Furthermore,  noise  in  tbe  measured  optical  flow  i_ 
accentuated  by  differentiation. 

In  tbe  continuous  approach,  the  whole  optical  flow 
field  is  used.  A major  shortcoming  of  both  the  local  and 
differential  approaches  is  that  neither  allows  for  errors  in 
the  optical  flow  d?ta.  This  is  why  we  choose  the  continuous 
approach  and  devise  a least-squares  technique  to  determine 
tbe  motion  of  tbe  camera  from  tbe  measured  optical  flow. 
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(1) 


The  proposed  algorithm  takes  the  abundance  of  available 
data  into  account  and  is  robust  enough  to  allow  numerical 
implementation. 

Independently,  Pratdny  chose  in  (13)  a similar  ap- 
proach to  outs,  lie  also  proposes  the  use  of  a least-squares 
method  to  determine  the  motion  parameters  but  never  dis- 
cusses how  exactly  this  is  to  lie  done.  Consequently,  he 
does  not  show  whether  his  scheme  can  be  used  to  uniquely 
determine  the  motion. 


9z=-T-qxi. 

W>i  define  the  components  of  t and  (3  u: 

T~(U,V,W)T  a^[A,H,C)T.  (2) 

Thus  we  can  rewrite  (I)  in  component  form: 


2.  Technical  Prerequisites 

In  this  section  we  review  the  equations  describing  the 
relation  between  the  motion  of  a camera  and  the  optical 
Bow  generated.  We  use  essentially  the  same  notation  as  [9]. 
A camera  is  assumed  to  move  through  a static  environment. 
l,ct  a coordinate  system  A *,Y,2  be  fixed  with  respect  to 
the  camera,  with  the  2-axis  pointing  along  the  optical  axis. 
If  w«  wish,  we  eau  Hunk  of  the  environment  as  moving  iu 
relation  to  this  coordinate  system.  Any  rigid  body  motion 
can  be  resolved  into  two  factors,  a translation  and  a rota- 
tion. We  will  denote  by  '/’  the  translational  component  of 
the  motion  of  the  camera  and  by  3 its  augular  velocity  (see 
also  Figure  1 which  is  redrawn  from  [9]}.  Let  the  initan- 
taueous  coordinates  of  a point  P in  the  environment  be 
(-Y  ,Y,2). 


Figure  I.  Coordinate  Systems 


(Note  that  Z > 0 for  points  in  front  of  the  imaging  sys- 
- tem.)  Let  f be  the  vector  (X,  Y,  Z)r,  where  T denotes  the 
transpose  of  a vector,  then  the  velocity  of  P with  respect 
to  the  A',  Y,  Z coordinate  system,  is: 


x :•  **-u-dz+cy 

Y>  cn-V-CX  + AZ  (3) 

Z'm-W-AY  + BX. 


where 1 denotes  differentiation  with  respect  to  time. 

The  optical  flow  at  each  point  in  the  image  plane  it 
the  instantaneous  velocity  of  the  brightness  pattern  at  that 
point.  Let  (7,  y]  denote  the  coordinates  of  a point  in  the 
image  plane  (see  Figure  1).  Since  we  assume  perspective 
projection  hetween  an  object  point  P and  the  corresponding 
image  point  p,  the  coordinates  of  p are: 
x Y 

*“?  *~z-  M 

The  optical  flow,  denoted  by  (u,  v),  at  a (mint  (x,y)  is: 

unri'  V ■*  y\  (5) 


Differentiating  (4)  with  respect  to  time  and  using  (3)  we 
obtain  the  following  equations  for  the  optical  flow: 


u 


v 


X'  XZ' 

Z 2l 

(-§  - » + Cy)  -*(-y-  Ay  + fix) 

Y'  ' YZ' 

Z 2* 

(-|  - Cx  + A)  — y(-y  - Ay  + fix). 


(6) 


We  can  write  these  equations  in  the  form: 


u *B  U,  -f-  Ur  V a=  V,  -)-  V,  (7) 


where  (ut,v,)  denotes  the  translational  component  of  the 
optical  flow  and  (ur,i/r)  the  rotational  component: 

u,  — U^zW  U,  sa  Azy  — /f(x*  -j- 1)  Cy, 
vt  — ~V  "t yW  v,  = A{y 1 + I)  - Dxy  — Cx. 

(») 

So  far  we  have  considered  a single  point  P.  To  define 
the  optical  flow  globally  we  assume  that  P lies  oo  a surface 
defined  by  a function  2 = Z(X,Y)  wh:  ‘h  is  positive  for  all 
values  of  X and  Y.  With  v:;  surface  and  any  motion  of 
a camera  we  can  therefore  associate  a certain  optical  flow 
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Mill  we  My  that  tbe  surface  ami  the  motion  generate  Util 
optical  flow. 

Optical  flow,  therefore,  depends  upon  the  six  parameters 
of  motion  of  the  ct-mcra  anti  upon  the  surface  whose  images 
are  analysed.  Can  ell  these  unknowns  b*  uniquely  recap* 
tured  solely  from  optical  flow?  Tbe  answer  is  no.  To  see 
this,  consider  a surface  £3  which  is  a dilation  by  a factor  k 
of  a surface  Si.  Further,  let  two  motions  denoted  by  A/| 
and  Af,  have  the  same  rotational  component  and  let  their 
translational  components  be  proportional  to  each  other  by 
the  same  factor  k (we  will  say  that  Af|  and  Af,  are  timilar). 
Then  the  optical  flow  generated  by  S\  and  A/|  it  the  tame 
as  the  optical  flow  generated  by  Sj  and  A/a.  This  follows 
directly  from  the  definition  of  optical  flow  (8).  It  is  still  an 
open  question  whether  there  are  any  other  pairs  of  distinct 
surfaces  and  motions  which  generate  the  same  optical  flow. 

Determining  the  motion  of  a camera  from  optical  flow 
is  much  easier  if  we  arc  told  that  the  motion  is  purely 
translational  or  purely  rotational.  In  the  next  two  scctious 
we  will  deal  with  these  two  special  cases.  Then  we  shall 
analyse  the  com  where  no  a priori  assumptions  about  the 
motion  of  the  camera  are  made. 

3.  Translational  Casa 


or: 

{/,  V,  - iV,Wi  - yUi  IV,  + xyWi  W t 

= - xV\ IV,  - yLljWi  + xylV,tV,.  (14) 

Since  we  assumed  that  A\  and  7)  and  /f,  and  7,  generate 
the  same  optical  flow,  the  above  equation  must  hold  for  aii 
z and  y.  Therefore  the  following  equations  have  to  hold: 

UiVj  » t/,V i 

« -I'tlVt  (15) 

-UlWj  =. 

These  equations  can  be  rewritten  as: 

UixVi’.Wi  a Ui'.Vi'.Wi  (16) 

from  which  it  follows  that  Aj  is  a dilation  of  A\.  It  is  clear 
that  the  scaling  factor  between  A\  and  A*  (or  equivalently 
between  and  7,)  cannot  be  recovered  from  the  optical 
flow,  regardless  of  the  number  of  points  at  which  the  flow  is 
known.  By  uniquely  determining  the  motion  of  the  camera, 
wc  will  mean  that  the  motion  is  uniquely  determined  up  to 
a constant  scaling  factor. 

3.2.  Least-Squares  Formulation 


In  this  section  we  discuss  the  case  where  the  motion 
of  the  camera  <*  assumed  to  be  purely  translational.  As 
Itcforc,  let  f sa  ( U , V,  IV)  be  the  velocity  of  the  camera. 
Then  the  following  equations  hold  (sec  (8)): 

-l/+xW  -V  + yW 

u,  - v,  ^ • (#) 


3.1.  Similar  Surfaces  and  Similar  Motions 

It  will  be  sliowu  next  that  if  two  flows  gcurrate  tbe 
same  optical  flow,  and  wc  know  that  tbe  motions  arc  purely 
translational,  then  the  two  surfaces  arc  similar  and  tbe  two 
camera  motions  arc  similar.  Let  A\  and  A]  be  two  surfaces 
and  let  1\  » (f/» , V, , W, )T  aud  7,  « (U,,  V,,H',)r  dcCct 
two  different  motions  of  a camera,  such  that  A\  and  7\  and 
At  and  7,  generate  the  same  optical  flow: 


-1/,+zlV, 

-Vi+yWi 

U as  ■ — * 

Ax 

Zx 

-(/,  -f  sWj 

.,  .-Vz  + vW'z 

“ Zt 

Zj 

(10) 

(») 


Eliminating  Z\,  Aj,  u and  v from  these  equations  we  obtain: 

-Ui+*Wt  __  -Uj  + xWj  . . 

-Vi+yWi  " -Vt  + yWj'  K 1 

We  can  rewrite  this  equation  as: 

H>,+ zM'iX-Va  + yW,) 

« {-V7  + xVV,X-Vi+  VIVO,  (13) 


In  general,  the  direction  of  the  optical  flow  at  two 
points  in  the  image  plane  determine  the  inotiou  of  a camera 
in  pure  translation  uniquely.  There  is  a drawback  however 
to  utilising  so  little  of  the  available  information.  The  opti- 
cal flow  we  measure  is  corrupted  by  noise  and  it  is  desirable 
to  develop  a robust  method  which  takes  this  into  account. 
Thus  wc  suggest  using  a least-squares  method  (4],  (14]  to 
determine  the  movement  parameters  and  the  surface  (i.e., 
tbe  best  lit  with  respect  to  some  norm). 

For  the  following  we  assume  that  the  image  plane  is 
the  rectangle  zrj— u>,  t«»)  aud  yr|— h,h).  Tbe  same  method 
applies  if  the  image  has  some  <>*  her  shape.  (In  fact,  it  can 
be  used  on  sub-images  corresponding  to  individual  objects 
in  the  case  that  the  environment  contains  objects  which 
may  move  relative  to  one  another).  Furthermore  we  have  to 
assume  that  \/Z  is  a bounded  function  and  that  the  set  of 
points  where  I /A  is  discontinuous  is  of  measure  zero.  This 
condition  on  1 /A  assures  us  that  all  necessary  integrations 
can  be  carried  out.  We  wish  to  minimize  the  following 
expression: 


LSI. 


((«- 


=£±*r+(.- 


(It) 


A 


In  this  case  then,  wc  determine  the  best  Gt  with  respect  to 
tbe  MLj  norm  which  is  defined  as: 


II  /(x,y)  11=  / J W,y)?dxdy.  (18) 
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The  steps  in  the  least- squares  method  Me  as  follows:  First 
we  determine  'iiat  7,  which  minimises  the  integrand  of  (17) 
at  every  point  (x,  y).  Then  we  determine  the  value;  of  U,  V 
and  IV  which  minimise  the  integral  (17). 

Let  us  introduce  the  following  abbreviations: 

a=z-U  + xW  p **  —V  + ylV.  (ltf) 

Note  that  the  expected  Dow,  given  V,  V and  IV  is  simpij. 

fl  ra  ^ 0 “ f ’ (-10) 

Then  we  can  rewrite  (17)  as: 


LJ-W  K“-|),  + (v~|),|dldV-  (21) 


Flgurt  2.  Geometrical  Interpretation 


We  proceed  now  with  the  first  step  of  our  minimisation 
method.  Differentiating  the  integrand  of  (17)  with  respect 
to  Z and  setting  the  resulting  expression  equal  to  uro 
yields: 

(—f)£+(v-f)Jr-o.  («) 

Therefore  we  can  write  Z as: 


1 

UQ  -f  V/J ' 


(23) 


This  equation,  by  the  way,  imposes  a constraint  on  U,  V 
and  IV,  since  Z must  be  positive.  We  do  not  make  use  of 
V is  except  to  help  us  pick  amongst  two  opposite  solutions 
for  the  translational  velocity  later  on.  Note  that  now: 


u 


up  — va 
dTfpi 


up  — va 


(24) 


and  we  can  therefore  rewrite  (17)  as: 


(ud  — «<*)* 
a>+/3* 


dzdy. 


(25) 


It  should  be  clear,  by  the  way,  that  uniformly  scaling  U,  V 
and  IV  does  not  change  the  value  of  (25).  This  is  a reflection 
of  the  fact  that  we  can  determine  the  motion  parameters 
only  up  to  a constant  factor. 


Before  proceeding  with  the  second  step,  we  give  a 
geometrical  interpretation  in  Figure  2 of  what  we  have  so 
far.  Supple  that  the  motion  parameters  U,  V,  and  IV  are 
given.  At  any  given  point,  say  (z0,  y0),  optical  flow  depends 
not  only  upon  the  motion  parameters  but  also  upon  the 
value  of  Z at  that  point,  Zq  sey.  However,  the  direction 
of  (u,  v)  does  not  depend  upon  Z0.  The  point  (u,  u)  must 
lie  along  the  line  L in  the  uv*plane  defined  by  the  equation 
uP~vaa=  0.  Let  the  menured  optical  flow  at  (7g,y0)  be 
denoted  by  (um,vm),  and  let  the  closest  point  on  the  line 
L be  (u,,  v,).  This  corresponds  to  a particular  Z,  given 
by  (23).  The  remaining  error  is  the  distance  between  the 
point  (um,vm)  and  the  line  L.  The  squme  of  this  distance 
is  given  by  the  integrand  of  (25). 

For  the  second  step,  we  differentiate  (25)  with  respect 
t i U,  V and  IV  Mid  set  the  resulting  expressions  equal  to 
sero: 


/:/: 


P[up  — vci)(uri  -f  vft) 


dzdy  0 


(•*+/»»)* 

fr  Jtiy  _ 


-kJ-w  (a»  + /?»)* 

Let  us  introduce  the  following  abbreviation: 

,,  _ (up  — VOr)(utt  + vP) 
(<*’ + />*)’  ‘ 

Then  equations  (26)  can  be  rewritten  as: 


(27) 


f J ((-V  + yW)K)dxdy  = 0 

-CL  ((—  V + xlV)ff ) dzdy  = 0 (28) 


i. 
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f f (I  - yV  I*  zV)h’\dxdy  0. 

The  sum  of  U tune*  the  first  integral,  t'  tunes  the  second 
integral,  ami  IF  tunes  (he  ttnnl  mt>  «ral  is  uh’iitiMlIy  lero. 
Thus  (he  three  equations  are  linear^  dependent.  j his  is  to 
lx  expected,  for  if: 


f[kU,kV,k\V)  »•  f(U,VtW), 


(2t>) 


where  / is  a dilfercntialde  function  and  k a constant,  then: 

0.  (30) 


U^J.  q.  vlL  .p  \y2L 

0U  dV  r 0W 


The  result  is  also  consistent  with  the  fact  that  only  two 
equations  arc  needed,  since  the  translational  velocity  can 
be  determined  only  up  to  a constant  factor.  Unfortunately 
equations  (28)  are  nonlinear  in  U,  V*  and  It'  and  we  are  not 
able  to  show  that  they  have  a unique  (up  to  a constant 
scaliug  factor)  solution. 

3.3.  Using  a Dillorent  Norm 

There  is  a way,  however,  to  devise  a least-squares 
method  which  allows  us  to  display  a closed  form  solution 
for  the  motion  parameters.  Instead  of  iiiiiimiizitig  (17),  we 
will  try  to  minimize  the  following  expression: 


/:/>- 


-u  + xty,, 


H-(o~ 


-V’l-yWn 


% 


)’l 


X («7  /?a)  ./x  Jy  (31) 

obtained  by  multiplying  the  integrand  of  (17)  by  a1  • 

Then  we  apply  the  same  least-squares  method  as  before 
to  (31).  When  the  measured  optical  flow  is  not  corrupted 
by  noise,  both  (31)  and  (17)  can  be  made  equal  to  icro  by 
substituting  the  correct  motion  parameters.  We  thus  obtain 
the  same  solution  for  the  motion  parameters  whether  we 
minimize  (31)  or  (17).  If  the  measured  optical  flow  is  not 
exact,  then  using  expression  (31)  for  our  minimization,  we 
obtaiu  the  best  fit  with  respect  not  to  the  M!  ? norm,  but 
to  another  norm  which  we  call  the  M Lap  norm: 

II  /(X, y}  ||„a=  f J |/(x,  y))*(ol  + P7) dx  dy.  (32) 

What  we  have  here  is  a minimisation  in  which  the  error 
contributions  are  weighted,  greater  importance  being  given 
to  points  where  the  optical  flow  velocity  is  larger.  This  is 
most  appropriate  when  the  measurement  of  larger  velocities 
it  more  accurate. 

Which  norm  gives  the  best  result*  depends  on  the 
properties  of  the  noise  in  the  measured  optical  flow.  The 
first  norm  is  better  suited  to  tbe  sitation  where  the  noise  in 
the  measurements  is  independent  of  the  magnitude  of  the 
optica)  flow.  Note  also  that  if  we  really  want  the  minimum 
with  respect  to  the  MLj  norm,  we  can  use  the  results  of  the 
minimization  with  respect  to  the  MLap  norm  as  stsrting 


values  hi  a numerical  minimization. 

We  discuss  now  our  least-squares  method  m the  case 
where  the  norm  is  chosen  to  lx  First  we  determine 

?.  by  differentiating  the  iutrgrnml  of  (31)  with  respect  to  Z 
and  setting  the  result  equal  to  zero.  Wc  again  get  (22): 


(u  _ ub. 


from  which  it  follows  that  (23): 


P, £ 

Vzi 


0, 


(33) 


*«£ i±£. 

un  -|-  up 


So  we  want  to  mutiniuc: 


f f (up  — vn)7  dxdy. 
J —kJ — *» 


(31) 


(35) 


(36) 


(38) 


Let  us  cat)  this  integral  y(U,  U,  IF),  then,  since: 
up  — uft  — ( vU  — uV)  — (xv  — yu)tt', 

we  have: 

9(U,V,W) 

« al/7  4-  M'*  -F e)V7  -i-  2 dUV  + 2«VIV  -)■  2JWU,  (37) 
where:  ^ 

on  f u1  dxdy 

fK  /“ 

4 = J J M*  dx  dy 
-IX  (xv  — yu)1  dx  dy 

d = — / / uvdxdy 

e~  / / u(xv  — yu) dxdy 

f — — J j v(xv  — yu)  dx  dy. 

Now  y(U,V,  W)  cannot  be  negative,  and  g(f/,  V,  IV)  = 0 
for  V — V ==  W = 0.  Thus  a minimum  can  be  found 
by  inspection,  but  is  not  what  we  might  have  hoped  for.  In 
fact,  to  determine  the  translational  velocity  using  our  least- 
squares  method  we  have  to  solve  the  following  homogeneous 
equation  for  t: 

Gf  = 0 (39) 

where  G is  tbe  matrix: 

fa  d A 

C = Jd  4 ej.  (40) 

Clearly  (39)  has  a solution  other  then  zero  if  and  only  if 
the  determinant  of  G is  zero.  Then  tbe  three  equations 
(39)  are  linearly  dependent  and  f can  be  determined  up 
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lo  * constant  factor.  In  general,  however,  as  the  data  is 
corrupted  by  noiiic,  9 cannot  be  made  equal  to  mo  for 
nonzero  translational  velocity  and  tot  ^ (0,0, 0)T  will  be 
the  only  solution  to  (30).  To  see  this  in  another  way,  note 
that  v has  the  following  form: 

g[W,kY,kW)  . k7g(U,V,W)  (II) 

where  k l*  a constant.  Clearly  F,  It')  assumes  its  min* 
itnal  value  for  U - V It'  — : 0. 

What  we  arc  really  interested  in,  is  determining  the 
direction  of  t which  minimises  9,  for  a fixed  length  of  t. 
lienee  we  impose  the  constraint  that  t lie  a unit  vector. 
If  t is  constrained  to  have  unit  magnitude,  the  minimum 
value  of  g is  the  smallest  eigenvalue  of  the  matrix  C and 
the  value  of  f for  which  9 assumes  its  minimum  can  he 
found  by  determining  the  eigenvector  corresponding  to  this 
eigenvalue  (tt).  This  follows  from  the  observation  that  9 is 
a quadratic  form  which  can  be  written  as: 

g(V,V,W)~tTGT.  (42) 

Note  that  G is  a positive  semidefinite  hcrmilian  matrix 
as  a 2!  0,  b j>  0,  c > 0,  oh  £ d*.  be  > e*  and  c«  > 
/*.  (The  last  three  inequalities  follow  from  the  Cauchy* 
Schwars  inequality  (T|,  [«]),  Hence  all  eigenvalues  are  real 
and  non*ncgative  ami  arc  the  solutions  X of  the  third  degree 
polynomial: 

X1 

— (a  -\-  b f c)XJ 

*4*  (n6  4*  be  4*  co  — d7  — • e*  — /*)X 

4-  (ae1  + bf7  + c d7  —abc  — 2 dtf)  0. 

There  is  an  explicit  formula  for  the  least  positive  root  in 
tcims  of  the  real  and  imaginary  parts  of  the  roots  of  the 
quadratic  resolvent  of  the  cubic.  In  our  ease  this  gives  us 
the  desired  smallest  root,  since  the  roots  cannot  be  negative. 
For  tbe  sake  of  completeness,  however,  various  pathologi* 
cal  cases  that  might  come  up  will  be  discussed  next,  even 
though  they  are  of  little  practical  interest. 

Note  that  X *=  0 is  an  eigenvalue  if  and  only  if  C is 
singular,  that  is,  if  the  constant  term  in  the  polynomial 
(•13)  equals  sero.  In  fact,  if  the  determinant  of  G is  scro 
one  can  Gnd  a velocity  t which  makes  9 sero.  It  follows 
from  a theorem  in  calculus  that  this  happens  only  when  the 
optical  How  is  either  not  corrupted  by  noise  at  all  or  only 
at  a few  points.  The  theorem  states  that  if  the  integral  of 
the  square  of  a bounded  and  continuous  function  is  sero 
then  the  function  itself  is  sero.  Hence  errors  can  only  occur 
at  points  where  the  optical  (low  is  discontinuous,  and  these 
are  exactly  the  points  where  the  surface  defined  by  Z is 
discontinuous.  (These  arc  also  the  places  where  existing 
methods  for  computing  the  optical  Sow  [5]  arc  subject  to 
large  errors). 

It  is  impossible  for  exactly  two  eigenvalues  to  be  sero, 
• this  would  imply  that  the  coefficient  of  X in  the  poly* 


nomi.il  (13)  equalled  sero.  while  the  cocllicicnt  of  X5  did 
not.  That  in  turn  would  imply  that  ab  =r  d7,  be  = e7,  and 
cn  f7,  while  u,  t,  ami  c arc  not  all  sero.  For  cqualit"  to 
hold  in  the  Cauchy-Scltvau  inequalities,  however,  u and  v 
must  Iwlli  l*e  proportional  to  xu  — yu.  This  can  only  be 
ttuc  (for  all  x ami  y in  the  image)  if  u«u^:0.  Hut  then  all 
six  integrals  become  sero  ami  consequently  all  three  eigen* 
values  arc  sero.  This  situation  is  of  little  interest,  since  it 
occurs  only  when  the  optical  How  data  is  sero  everywhere. 
Then  (he  velocity  is  zero  too.  Once  the  smallest  eigen- 
value is  known,  it  is  straightforward  to  find  the  translational 
velocity  which  best  matches  the  given  data.  To  determine 
the  eigenvector  corresponding  to  an  eigenvalue,  say  X],  wt 
have  lo  solve  the  followiug  set  of  liuear  equations: 

(a  — X,)tf  4*  dV  + JW=*  0 

dU  + (b-\i)V+tW**0  (44) 

JU  + cF4-(c-X,)IV-0. 

As  Xt  is  an  eigenvalue,  equations  (44)  arc  linearly  depen- 
dent. Let  us  for  a moment  assume  that  all  eigenvalues  we 
distinct,  that  if,  tbe  rank  of  the  matrix  (G  — X/),  where  / 
is  the  identity  matrix,  is  two.  Then  we  can  use  any  pair  of 
them  to  solve  for  U,  V in  terms  of  IF  say.  There  are  three 
ways  of  doing  this.  For  numerical  accuracy  we  may  add  the 
three  results  to  get  the  symmetrical  forms: 

(/  aa  (h  — X|Xe  X|)  — /(b  — Xi) 

— «f(c  — X 1 ) I*  t(f  4-  d — e) 

F (c  — X,)(n  — Xj)  — d[c  — Xt)  . . 

— 1(«  — X|)  4-  /(</  4*  < — /) 

IF  ta  (a  — X|)(i  — X|)  — e(«  — Xj) 

-/(b-X,)4-d(e4-/-<0* 

Note  that  X|  will  be  very  small,  if  the  data  is  goad,  and 
one  may  wish  to  just  approximate  the  exact  solution  by 
using  the  nlmve  equations  with  Xi  set  to  scro.  (Then  there 
is  no  need  to  find  the  eigenvalue).  In  any  case,  the  result- 
ing velocity  may  now  be  normalised  so  that  its  magnitude 
equals  one.  There  is  one  remaining  dilliculty,  arising  from 
the  fact  that  if  t is  a solution  to  our  minimization  problem, 
so  is  —t.  Only  one  of  these  solutiou  will  correspond  to 
positive  values  of  Z in  equation  (31)  however.  This  can  be 
ensily  seen  by  evaluating  (31)  at  some  point  in  the  image. 
The  ease  where  the  two  smallest  eigenvalues  arc  the  same 
will  be  discussed  in  one  of  the  uexl  paragraphs. 

There  is  a simple  geometrical  interpretation  of  what 
we  have  done  so  far.  To  this  end  we  consider  the  surface 
defined  by  g{U,  F,  IF)  — k where  it  is  a constant.  Note 
that  we  can  always  Gnd  a new  coordinate  system  U,V,W 
in  which  g[U,  V,  W)  can  be  written  as: 

X|lf*  4*  X|F*  4-  XjlF*  = k (46) 

where  X,  for  t = 1,2,3  are  tbe  three  eigenvalues  of  the 
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quadratic  form.  If  the  eigenvalues  arc  all  non-zero,  the 
surface  g{U,  V,  M'}  *=  k is  au  ellipsoid  with  three  orthogonal 
semi-axes  of  length  y/k/X ,\  We  arc  particularly  interested 
in  the  case  where  the  constant  k is  the  smallest  eigenvalue. 
Then  all  three  sewi-aues  hare  lengths  less  than  or  equal 
to  one.  lienee  the  ellipsoid  lies  within  the  unit  sphere.  If 
the  two  smallest  eigenvalues  arc  distinct,  the  unit  sphere 
touches  the  ellipsoid  in  two  places,  corresponding  to  the 
largest  axis.  If  the  two  smaller  eigenvalues  happen  to  be  the 
same,  however,  the  unit  sphere  touches  the  ellipsoid  along 
a circle  and  as  a result  ail  the  velocity  vectors  lying  in  a 
plane  spanned  by  two  eigenvectors  give  equally  low  errors. 
Finally,  if  all  three  eigenvalues  arc  equal,  no  direction  for 
f is  preferred,  since  the  ellipsoid  becomes  the  unit  sphere. 

The  case  where  exactly  one  eigenvalue  is  zero  also  has 
a simple  geometrical  interpretation.  The  surface  defined  by 
9{U,\',W)  - 0 is  a straight  line,  which  can  Iw  seen  easily 
?;ont  the  following  equation: 

X»&*  + XiP*«0  (47) 

written  for  the  case  when  Xj  is  xero.  (Remember  that  X| 
and  Xj  arc  both  positive.)  Clearly  the  unit  sphere  intersects 
this  line  in  exactly  two  points,  one  of  them  corresponding 
to  imsitivc  values  for  V,  in  equation  (34). 

The  method  which  we  just  described  can  be  easily 
implemented.  To  this  cud,  the  problem  can  be  discrctixed. 
An  expression  similar  to  (31)  can  be  derived  where  the 
integrals  arc  approximated  by  sums.  Our  minimisation 
method  cau  then  tic  applied  to  these  sums.  The  result- 
ing equations  are  similar  to  ones  descrilicd  m this  section, 
wilh  summation  replacing  integration.  We  implemented 
the  resulting  algorithm  and  tested  it  using  synthetic  data 
including  additive  noise.  The  results  agreed  with  our  ex- 
pect ations. 

One  can  use  the  ratio  of  the  l»sgcst  to  the  small- 
est eigenvalue,  the  so  called  condition  number  (15),  as  a 
measure  of  confidence  in  the  computed  velocity.  The  result 
it  very  sensitive  to  errors  iu  the  measurements  unless  this 
ratio  is  much  bigger  than  one. 

Curiously,  the  same  error  integral  at  (35)  it  obtained 
in  the  case  where  the  Af  Lz*v  norm  is  used: 

II  /(*. »)  ||zwv  /_4/_  J/(I.V)Z(*,V)),(«1  + v*)dxdy. 

(<8) 

We  can  arrive  at  a similar  solution  by  multiplying  the  in- 
tegrand in  (17)  by  X7  instead  of  by  a7  + 07.  In  that  case 
the  minimisation  is  carried  out  with  respect  to  the  MLx 
norm  defined  by: 

II  /(*.  v)  ||z=  / f |/(x, y)X(x, y))1  dx  dy.  (49) 


We  end  up  with  a quadratic  form  similar  to  g,  only  the 
integrals  for  the  six  constants  corresponding  to  a,  b,  c,  d, 
e,  and  / arc  a bit  more  complicated.  Curiously  they  only 
depend  on  the  direction  of  the  optical  How  at  each  point, 
not  its  magnitude. 

Also,  other  constraints  could  he  used.  If  we  insist 
on  U7  -|-  V’  = I,  for  example,  we  obtain  a quadratic 
instead  of  a cubic  equation,  and  if  we  use  H'  1,  a linear 
equation  only  need  to  lie  solved.  The  disadvantage  of  these 
approaches  is  that  the  result  is  sensitive  to  the  orientation  of 
the  coordinate  axes.  Clearly,  in  the  case  of  exact  data,  we 
get  the  right  solution  using  any  of  the  constraints  mentioned 
above. 


3.4.  Using  a Different  Constraint 

The  minimization  scheme  discussed  in  the  previous  sec- 
tion gives  us  a unique  solution  in  most  cases  for  the  velocity 
vector  t.  Here  we  promise  a slightly  different  approach 
which  always  gives  us  a unique  solution.  Note  that  apply- 
ing the  first  step  m our  iiii'..ittii>alion  method  gives  us  a 
'eonitraiut  between  the  values  of  a,  the  velocity  vector  and 
tlie  fijilical  flow  at  evciy  |>oiul.  We  can  in  addition  assume 
that  2 a*  7.o  is  known  at  a particular  point,  say  (xg,yo). 
Using \tlie  A //»s„v  norm  iu  our  scheme,  we  want  to  mini- 
mite: v 

fJjuX  - (-U  + *IV)|*  + \v7.  - {-V  + y »V)1* 

X (u1  + u*)dzdy.  (50) 

Differentiating  (50)  wilh  respect  to  %,  and  setting  the  result- 
ing expression  equal  to  tcro,  we  obtain: 


UOr  Vp 
U*  -f-  V*  ’ 


(51) 


Thus  we  propose  to  solve  the  minimisation  scheme  under 
the  following  constraint: 


20(u§  + vj)  — (uoa0  -f  vtPo)  0 (52) 


where  Ug  and  v0  denote  the  components  of  the  optical 
Bow  measured  at  (z0,yo)  and  a0  and  0O  denote  a and  fi 
evaluated  at  (xo,yo).  The  error  integral  (50)  becomes  after 
substituting  (51): 

/:/:  (u0  — vet)1  dx  dy  (53) 

which  is  the  same  as  (35)  and  is  denoted  by  g[U,  V,  W)  (37). 
Thus  we  want  to  ininimiie: 


y(l/ , V,  W)  -f  2X[#o(«/q  + Ug)  — (uo«g  + wo0o))  (54) 


Here  optical  flow  velocities  for  points  which  are  further 
away  are  weighted  more  heavily.  This  is  most  appropriate 
when  the  measurement  of  larger  velocities  is  Jess  accurate. 


where  X denotes  a Lagrangiau  multiplier.  To  determine 
U,V  aud  W the  following  linear  equations  obtained  by 
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differentiating  (51)  with  rcsjicct  to  U,V,W  ami  X have  to 
be  solved; 


alf  + dV  -|-  /IV  + Xu0  = 0 

dU  4-  6V  + cIV  *(•  Xi-o  = 0 

fU  cV  + cW  - X(x„«o  + yot’c) » 0 

uoU  + — (*o“o  + yoi*a)IV  *»  — Zo(uJ  + vl)‘ 

Thcic  equations  can  be  written  in  the  form: 

Cx?x  - P 


(M) 


(56) 


where  - (U,V,W,\)T  and  F = (0,0/d, -Z0(“o  + 
vJ))T.  Let  the  determinant  of  G*  be  Ao: 


A0  «■  (d*  - o6X*o«o  + Wo)* 

+(e,-hcK+(/1-ne)v| 
+3((de  - bf)  M0(*0«e  + Vo«o) 
+(d/  - ae)v0(*oUo  + KoVo) 
+(«f  — e/)«pVo). 


(57) 


Assuming  that  Ao  ft  0 we  can  easily  determine  T\  from 
(55): 

Tx-g^p.  («») 

Introducing  the  following  abbreviation: 

^o(«o  + vo) 


Ao 


(so) 


we  can  give  these  formulae  for  t\'. 

U ■*  A'(uo(hc  — e1)  -)•  «o(«/  — «<0 

+ (*o«o  + »«Vo  XV  ~ *)1 

V — K(u0(e/  — cd)  -f  u0(oe  — /*) 

+ (*o«o  + Vo«oXa«  “ «*/)l  (6°) 

W«K|uo(de-k/)  + vo(d/-oe) 

+ (*ouo  + VoVoXd1  — «5)1 
X - Aloe*  + edl  + V* -ohc- 2de/j. 

The  disadvantage  of  this  approach  is  that  the  result  depends 
upon  Ur  values  of  the  opUcal  flow  at  a single  point.  To 
circumvent  this  problem  we  propose  to  determine  average 
values  for  U,V  and  W in  the  following  manner.  First  note 
that  we  are  only  interested  in  the  ratios  of  U/W  and  V/W 
which  obviously  do  not  depend  upon  the  (unknown)  value 
for  Z,.  Equivalently  we  could  determine  thr*  value  for  K 
from  the  condition  that  t should  have  unit  length.  Hence 
we  can  determine  values  for  (/,  V and  IV  which  depend  only 
upon  the  values  of  the  optical  flow  at  a single  point  and  the 
coefficient  in  the  matrix  G.  If  wc  want  to  remove  the  de- 
pendence or  the  result  on  the  data  at  a single  point,  we  can 
simply  average  the  values  obtained  for  all  image  points. 

In  the  cue  where  the  data  is  exact,  i.e.,  where  the 


determinant  of  G,  denoted  by  det  G,  vanishes,  the  solution 
for  f is  the  same  one  as  obtained  using  no  constraint  in 
our  minimisations  scheme.  To  sec  this  just  observe  that 
in  that  case  X =»  0 as  X «»  —K'dctG).  In  the  case  where 
A0  =*  0,  equations  (55)  have  a solution  only  when  detG 
0.  We  do  not  have  to  be  concerned  with  the  case  where 
A0  **  0 but  detG  0 as  we  can  argue  that  equations 
(5f>)  always  have  to  have  a solution.  Note  that  our  method 
i?  based  on  the  condition  that  is  a certain  function  (51) 
of  l/,V,W.  Hence  (52)  cannot  impose  a constraint  which 
would  be  impossible  to  satisfy. 

The  methods  discussed  in  this  section  have  been  ap- 
plied to  noisy  synthetic  data  with  the  expected  results. 

4.  Rotational  Casa 

Sup;>osc  now  that  the  motion  of  the  camera  is  purely 
rotational.  In  order  to  determine  the  motion  from  optical 
flow  we  again  use  a least-squares  algorithm  with  the  M L* 
norm  described  in  the  previous  section.  Recall  that  in  this 
case  the  optical  flow  is  (see  (8)): 


V,  : 


Axy-flfx’  + lJ  + Cy 

A(y*  + 1)  - Dxy  - Gx. 


(61) 


We  wilt  show  now  in  an  analogous  fashion  to  section  3.1 
that  two  different  rotations,  say  Oi  » (Ai,B\,C\)t  and 
wj  cs  (dj,  Wj,Cj)T,  cannot  generate  the  same  optical  flow. 
Assuming  the  converse,  the  following  equations  have  to  hold 
for  all  values  of  x and  y: 


toy  — Hi(x*-f  1)  + Cty 

* 'toy  - Aj(xa  + 1)  + Ci y 
Aj(y*-f  l)  — flixy  — Gj* 

*•  Aj(y*  -(- 1)  — JJjxy  — Cjx 


(62) 


from  which  we  can  immediately  deduce  that  «*  Qt. 

In  general,  the  direction  of  the  optical  flow  at  two 
points  and  its  magnitude  at  ouc  point  determine  the  motion 
of  a camera  in  pure  rotation  uniquely.  We  choose  instead 
to  minimize  the  following  expression: 


CL  |(u  — uf  )*  + (v  — v,)*)  dx  dy. 


(63) 


As  the  motion  is  purely  rotational,  the  optical  flow  does 
not  depend  upon  the  distance  to  the  surface  and  therefore 
we  may  omit  the  first  step  in  our  method.  Thus  we  im- 
mediately differentiate  (63)  with  respect  to  A,  B and  G and 
set  the  resulting  expressions  equal  to  sero: 

U')*V  + («- VrXy*  + 1)1  dx  dy  = 0 
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IX  l(u  “ *rX**  4*  0 4*  (»  — vr)*y)  itdysa  0 (64) 


J_J_  Ku  — “r)v  “ O'  ~ ur)ij  dz  dy  **  0. 


Rewriting  the  above  equations  we  obtain: 

fj  tu*V  + V(V*  4-  01  dxdy 


“ / J I«rxy  4-  v, (y1  + l)]  dz dy 

/:/:  (nj*1  4*  1)4*  v*y]dzdy 
= /_ A/_w  M**  *H)  + "riy)  <1* 

“»  / / (ury  — v.x]  ax  <fy 

* —kJ  —‘I* 


and  expanding  these  equations  yields: 
aA  + Xll  + ]C~t 
2a  ■)•  5/J  4*  ic  » 7 
7a  4*  4- 


where: 


and: 


a « j_J_J*7y7  + (v*  4-  0*1  dx  dy 

f kf  ((** + l)s  + i*y*ld*dy 

***  LJ„Jxi  + yi\dxdv 

3 “ - f_J_JxV{**  + y*  + 2)]  dx  dy 

7— /_£***- 
* “ f_J_Ju*v + Wn*  4-  0)  * <*v 


- J“(x*  + 0 4*  wxy)  dx  dy 


it»  — /:/>  — uxjdxdy. 


(65) 


(66) 


(67) 


(68) 


Thus,  provided  the  matrix  A/  is  non-singular,  we  can  com- 
pute the  rotation  as  follows: 

s3»M"*,it.  (70) 


It  is  easy  to  see  that  the  matrix  Af  is  non-singular  in  the 
special  case  of  a rectangular  image  plane  since  then: 


. j a*1*  j 1 i 

9 *»  4wh(~-  4-  — -{- 1)  4*  — g — 

I j uw*  i 2»»*  i ,\  , -l>"5As 

5 =*  4wh(y  -i — ^ — h 1)  4-  — — 
•ln»h 


(70 


(u^f-A*) 

«l 

3 = a « 7 “ o. 


So  in  the  case  of  a rectangular  image  plane,  the  matrix  is 
diagonal,  which  makes  it  particularly  easy  to  compute  its 
inverse.  In  fact,  the  matrix  is  diagoual  if  the  image  plans 
is  symmetrical  about  the  x-axis  and  the  y-axis.  As  the 
extent  of  the  image  plane  decreases,  however,  the  matrix  M 
becomes  ill  conditioned.  That  is  inaccuracies  in  the  three 
integrals  (£,  7,  and  r7i)  computed  from  the  observed  law 
arc  greatly  magnified.  This  makes  sense  since  we  cannot 
expect  to  accurately  determine  the  component  of  rotation 
about  the  optical  axis  when  observations  arc  confined  to  a 
small  cone  of  directions  about  the  optical  axis. 

Again,  in  our  numerical  implementation  of  tne  algo- 
rithm the  integrals  in  (67)  can  be  approximated  by  sums. 
The  methods  discussed  in  this  section  have  been  applied  to 
noisy  synthetic  data  with  the  expected  results. 


5.  General  Motion 


We  would  like  now  to  apply  a least-squares  algorithm 
to  determine  the  motion  of  a camera  from  optical  low 
given  no  a priori  assumptions  about  the  motion.  It  is 
plain  that  a least-squares  method  is  easiest  to  use  when  the 
resulting  equations  nre  linear  in  all  the  motion  parameters. 
Unfortunately,  there  exists  no  norm  which  will  allow  us 
to  achieve  this  goal.  We  did  find  a norm,  however,  which 
resulted  in  equations  that  are  linear  in  some  of  the  un- 
knowns and  quadratic  in  the  others.  We  again  attack  the 
minimitation  problem  using  the  norm  under  the 

constraint  that  (/*  4*  V7  4-  W*  « 1.  The  ensueing  equa- 
tions nre  polynomials  in  the  unknowns  U,V,W,A,n  and 
C and  can  be  solved  by  a standard  iteration  method  like 
Newton's  method  or  Bairstow’s  method  |14)  or  by  an  inter- 
polation scheme  like  rtguU  Msi  (H).  The  expression  we 
wish  to  minimise  is: 


If  we  call  the  coefficient  matrix  in  (66)  M and  the  column 
vector  on  the  right-hand  side  II,  then  we  have: 

MQ  = «.  (60) 


/_k/j|u“(  |+**0]*4-I»-(^4-v,))*Mo*4-/»*)^^. 

(72) 

The  first  step  is  to  differentiate  (be  integrand  of  (72)  with 
respect  to  Z and  set  the  resulting  expression  equal  to  sero: 
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nl  *|*  / 9* 

(•<  — u,)a  4*(u  ~vf)/T 


W«  introduce  the  Langrangun  multiplier  X a*  before  and 
attempt  <o  minirniic: 


/A 

J_v  ((“  -«,)/»-  (u  - VrJoJa  rfx  dy 

+ X(UI  + V'J  + »'*-!).  (74) 

The  equations  we  have  to  lotve  to  determine  the  motion 
parameters  are  obtained  by  differentiation: 

X (-- xy/J  -l*  (y3  |*  l)n]dxdy  <-  0 

/ J-  J^U  “ u')^  ~ (u  ~ 

X ((x3  -|*  l)/f  — xyojdr  dy  s*  0 

j J _wl{u  ~ M')0  * (v  - »»W 

X [y/J  *!-  xnjj  dx  dy  * 0 

J J 

X (w  — of)  dx  «fy  •!*  \U  »«  0 

/.a/.w((u  “ Uf  ^ ~ (v  “ U'W 

X (««  — «,)dxdy  — XV  *»  0 

/JjN  ~ u-)^  - (®  - 

>•  ((“  - Wr)v  + (»  — Vr)xJ  dx  dy  XIV  w 0 

V3  + V1  *}•  IV3  mi, 

(75) 

Note  that  the  first  three  of  these  equations  arc  linear  in 
A I)  and  C from  which  these  parameters  can  lie  determined 
uniquely  in  terms  of  U,  V aud  IV.  Then  we  can  determine 
V,  V’  aud  IV  from  the  last  four  equations  by  a numerical 
method.  To  this  end,  the  problem  can  be  discretized  and 
equations  analogous  to  (75)  derived,  where  summation  of 
the  appropriate  expressions  is  used  instead  of  integration. 

0.  Summary 

Our  objective  was  to  devise  a method  for  determin- 
ing the  motion  of  a camera  from  optical  Bow  which  allows 
for  noise  in  the  measured  data.  The  least-squares  method 
which  we  proposed  in  this  paper  meets  our  goal  and  n 
also  suitable  for  numerical  implementation.  An  important 


application  of  our  results  is  in  passive  navigation.  Here 
the  path  and  instantaneous  altitude  of  a vehicle  it  to  be 
determined  from  information  gleaned  about  the  cnviion- 
ment  without  the  emission  of  sampling  radiation  from  the 
vehicle. 
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